clear all;
[A, B, E, H, N] = readD();
eqn.A_=A;
eqn.B=B;
eqn.C=B';
eqn.N_=N;
opts.blyap.maxiter=20;
% Set tolerances and iteration limit for the Lyapunov-plus-positive solver
opts.blyap.res_tol  =  1e-6;
opts.blyap.rel_diff_tol  =  1e-6;
opts.blyap.maxiter = 6;

% Set options for the ADI-based inner Lyapunov solver
% and related shift computation
opts.shifts.num_Ritz = 50;
opts.shifts.num_hRitz = 25;
opts.shifts.method = 'projection';
opts.shifts.num_desired = 6;

opts.adi.maxiter = 24;
opts.adi.res_tol = 1e-12;
opts.adi.rel_diff_tol = 0;
opts.adi.info = 2;
opts.adi.norm = 'fro';
opts.norm = 'fro';

% Set mess_res2_norm options
opts.resopts.res.maxiter = 8;
opts.resopts.res.tol = 1e-6;
opts.resopts.res.orth = false;

opts.fopts.LDL_T = false;
opts.fopts.norm = 'fro';

opts.srm.tol = 1e-10;
opts.srm.info = 1;
[oper, opts] = operatormanager(opts, 'default');
[out,~,~,~]=mess_lyapunov_bilinear(eqn,opts,oper);
%[out,~,~,~]=mess_lyapunov_bilinear(eqn);
V1 = out.Z;
V = V1*V1';